{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3261888c",
   "metadata": {},
   "source": [
    "### Selection of misexpression-associated SVs \n",
    "\n",
    "1. **Selection of misexpression-associated and control SVs**\n",
    "    * Misexpression variants: \n",
    "        * To enrich for variants associated with misexpression, we selected variants with a median TPM > 0.5 and median z-score > 2 while removing variants where a carrier has expression < 0.1 TPM. \n",
    "        * This retains variants with some variability around the misexpression cutoff but removes variants with bimodal expression distributions that are likely non-causal. \n",
    "\n",
    "    * Control variants:\n",
    "        * Selected variants in +\\- 200kb in cis to gene that were misexpressed across all INTERVAL samples\n",
    "        * Selected variants with a maximum TPM = 0 across all inactive genes within 200kb\n",
    "    \n",
    "2. **Annotation of misexpression-associated SVs with VEP**\n",
    "3. **Annotation of misexpression-associated SVs by position**\n",
    "4. **Comparison of misexpression-associated and control SV length**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "26432c81",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys \n",
    "import pandas as pd \n",
    "import pysam\n",
    "import numpy as np\n",
    "from pathlib import Path\n",
    "from collections import Counter\n",
    "from scipy.stats import mannwhitneyu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2ae029b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "wkdir=\"/lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/misexpression_v3/\"\n",
    "wkdir_path = Path(wkdir)\n",
    "# inputs \n",
    "sv_info_path = \"/lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/lof_missense/data/sv_vcf/info_table/final_sites_critical_info_allele.txt\"\n",
    "express_info_dir = wkdir_path.joinpath(\"4_vrnt_enrich/sv_count_carriers/gene_body/200kb_window/express_carrier_info\")\n",
    "misexp_genes_path = wkdir_path.joinpath(\"2_misexp_qc/misexp_metrics/misexp_genes_tpm0.5_z2.txt\")\n",
    "vep_msc_path = wkdir_path.joinpath(\"4_vrnt_enrich/sv_vep/msc/SV_vep_hg38_msc_parsed.tsv\")\n",
    "vep_all_path = wkdir_path.joinpath(\"4_vrnt_enrich/sv_vep/all/SV_vep_hg38_all_parsed.tsv\")\n",
    "gencode_path = wkdir_path.joinpath(\"reference/gencode/gencode.v31.annotation.sorted.gtf.gz\")\n",
    "# output \n",
    "output_dir = wkdir_path.joinpath(\"5_misexp_vrnts/test_cntrl_sets/\")\n",
    "output_dir.mkdir(parents=True, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "12f3e090",
   "metadata": {},
   "outputs": [],
   "source": [
    "# constants \n",
    "SV_TYPES = [\"DEL\", \"DUP\", \"INV\", \"MEI\"]\n",
    "CHROMOSOMES = ['chr1', 'chr2', 'chr3', 'chr4', 'chr5', 'chr6',\n",
    "               'chr7', 'chr8', 'chr9', 'chr10', 'chr11', 'chr12',\n",
    "               'chr13', 'chr14', 'chr15', 'chr16', 'chr17', 'chr18',\n",
    "               'chr19', 'chr20', 'chr21', 'chr22']\n",
    "# MAF cutoffs \n",
    "af_lower = 0\n",
    "af_upper = 0.01\n",
    "# misexpression cutoffs \n",
    "z_cutoff_misexp = 2 \n",
    "tpm_cutoff_misexp = 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "52d68f9b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of misexpressed genes: 4437\n"
     ]
    }
   ],
   "source": [
    "# load structural variant info\n",
    "sv_info_df = pd.read_csv(sv_info_path, sep=\"\\t\", dtype={\"plinkID\": str}).rename(columns={\"plinkID\":\"vrnt_id\"})\n",
    "sv_info_id_af_df =  sv_info_df[[\"vrnt_id\", \"AF\"]]\n",
    "# misexpressed genes\n",
    "misexp_genes = pd.read_csv(misexp_genes_path, header=None)[0].unique()\n",
    "print(f\"Number of misexpressed genes: {len(misexp_genes)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b238b658",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chr1\n",
      "chr2\n",
      "chr3\n",
      "chr4\n",
      "chr5\n",
      "chr6\n",
      "chr7\n",
      "chr8\n",
      "chr9\n",
      "chr10\n",
      "chr11\n",
      "chr12\n",
      "chr13\n",
      "chr14\n",
      "chr15\n",
      "chr16\n",
      "chr17\n",
      "chr18\n",
      "chr19\n",
      "chr20\n",
      "chr21\n",
      "chr22\n"
     ]
    }
   ],
   "source": [
    "# combine all variant-gene-sample pairs in 200kb windows \n",
    "all_sv_in_window_df_list = []\n",
    "misexp_sv_gene_list, cntrl_sv_gene_list = [], []\n",
    "express_info_path = Path(express_info_dir)\n",
    "for chrom in CHROMOSOMES:     \n",
    "    print(chrom)\n",
    "    sv_intersect_express_path = express_info_path.joinpath(f\"{chrom}_express_carrier_info.tsv\")\n",
    "    sv_intersect_express_df = pd.read_csv(sv_intersect_express_path, sep=\"\\t\")\n",
    "    sv_intersect_express_genes_df = sv_intersect_express_df[sv_intersect_express_df.gene_id.isin(misexp_genes)]\n",
    "    # add SV information \n",
    "    sv_intersect_express_info_df = pd.merge(sv_intersect_express_genes_df, \n",
    "                                            sv_info_id_af_df,\n",
    "                                            how=\"inner\", \n",
    "                                            on=[\"vrnt_id\", \"AF\"])\n",
    "    sv_in_window_df = sv_intersect_express_info_df[(sv_intersect_express_info_df[\"AF\"] >= af_lower) &\n",
    "                                                   (sv_intersect_express_info_df[\"AF\"] < af_upper) & \n",
    "                                                   (sv_intersect_express_info_df.genotype.isin(['(0, 1)', '(1, 1)']))].copy()\n",
    "    # add variant-gene ID  \n",
    "    sv_in_window_df[\"vrnt_gene_id\"] = sv_in_window_df.vrnt_id + \"|\" + sv_in_window_df.gene_id \n",
    "    # add rare SVs in cis to misexpressed genes \n",
    "    all_sv_in_window_df_list.append(sv_in_window_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c9c293c6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of unique variants in windows: 23159\n",
      "Number of control variants: 20150\n",
      "Max TPM of gene in cis to controls: 0.0\n",
      "Max z-score of gene in cis to controls: -0.0187738852838987\n"
     ]
    }
   ],
   "source": [
    "# combine all chromosomes  \n",
    "all_sv_in_window_df = pd.concat(all_sv_in_window_df_list)\n",
    "# all variants in 200 kb windows\n",
    "vrnt_id_in_windows = set(all_sv_in_window_df.vrnt_id.unique())\n",
    "print(f\"Total number of unique variants in windows: {len(vrnt_id_in_windows)}\")\n",
    "\n",
    "# control variants: all variants with max expression of gene in cis window is TPM = 0 \n",
    "all_sv_in_window_max_tpm_df = all_sv_in_window_df.groupby(\"vrnt_id\", as_index=False).TPM.max()\n",
    "vrnt_gene_id_cntrl_df = all_sv_in_window_max_tpm_df[all_sv_in_window_max_tpm_df.TPM == 0]\n",
    "vrnt_id_cntrl = set(vrnt_gene_id_cntrl_df.vrnt_id.unique())\n",
    "print(f\"Number of control variants: {len(vrnt_id_cntrl)}\")\n",
    "# checks\n",
    "print(f\"Max TPM of gene in cis to controls: {all_sv_in_window_df[all_sv_in_window_df.vrnt_id.isin(vrnt_id_cntrl)].TPM.max()}\")\n",
    "print(f\"Max z-score of gene in cis to controls: {all_sv_in_window_df[all_sv_in_window_df.vrnt_id.isin(vrnt_id_cntrl)]['z-score'].max()}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b6ffefcf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of variants passing quantile cutoff: 124\n",
      "Number of misexpression variants: 105\n",
      "Overlap between control and misexpression SVs: 0\n"
     ]
    }
   ],
   "source": [
    "# misexpression SVs, median TPM > 0.5 and z-score > 2\n",
    "vrnt_gene_median_tpm_z_df = all_sv_in_window_df.groupby([\"vrnt_gene_id\"], as_index=False)[[\"TPM\", \"z-score\"]].quantile(q=0.5, interpolation='linear')\n",
    "vrnt_gene_median_tpm05_z2 = set(vrnt_gene_median_tpm_z_df[(vrnt_gene_median_tpm_z_df.TPM > tpm_cutoff_misexp) & \n",
    "                                                          (vrnt_gene_median_tpm_z_df[\"z-score\"] > z_cutoff_misexp)].vrnt_gene_id.unique())\n",
    "vrnt_id_pass_median = set([vrnt_gene_id.split(\"|\")[0] for vrnt_gene_id in vrnt_gene_median_tpm05_z2])\n",
    "print(f\"Number of variants passing quantile cutoff: {len(vrnt_id_pass_median)}\")\n",
    "# remove SVs with a carrier with TPM < 0.1 \n",
    "vrnt_gene_id_median_to_rmv = set(all_sv_in_window_df[(all_sv_in_window_df.vrnt_gene_id.isin(vrnt_gene_median_tpm05_z2)) & \n",
    "                                                     (all_sv_in_window_df.TPM < 0.1)].vrnt_gene_id.unique())\n",
    "vrnt_gene_id_median_pass = vrnt_gene_median_tpm05_z2 - vrnt_gene_id_median_to_rmv\n",
    "vrnt_id_misexp = set([vrnt_gene_id.split(\"|\")[0] for vrnt_gene_id in vrnt_gene_id_median_pass])\n",
    "print(f\"Number of misexpression variants: {len(vrnt_id_misexp)}\")\n",
    "# check overlap between control and misexpression variants \n",
    "print(f\"Overlap between control and misexpression SVs: {len(vrnt_id_cntrl.intersection(vrnt_id_misexp))}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "9c22a0ed",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of misexpression variant carriers: 150\n"
     ]
    }
   ],
   "source": [
    "# total variant carriers pairs\n",
    "misexp_vrnt_carrier_df = all_sv_in_window_df[all_sv_in_window_df.vrnt_gene_id.isin(vrnt_gene_id_median_pass)]\n",
    "print(f\"Number of misexpression variant carriers: {len(misexp_vrnt_carrier_df)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e4fb4aaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# write variant carrier information to file \n",
    "misexp_vrnt_carrier_info_df = misexp_vrnt_carrier_df[[\"vrnt_id\", \"egan_id\", \"rna_id\", \"genotype\", \"gene_id\", \"TPM\", \"z-score\"]]\n",
    "misexp_vrnt_carrier_info_path = output_dir.joinpath(\"misexp_vrnt_gene_smpl.tsv\")\n",
    "misexp_vrnt_carrier_info_df.to_csv(misexp_vrnt_carrier_info_path, sep=\"\\t\", index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a8cdf64a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# misexpression-associated gene-SV pairs \n",
    "misexp_sv_gene_pairs_df = misexp_vrnt_carrier_df[[\"vrnt_id\", \"gene_id\"]].drop_duplicates()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "292839c4",
   "metadata": {},
   "source": [
    "**Count number of SV types in control and misexpression groups**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "9d8fba94",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of SVs in controls and misexpression groups: 20255\n"
     ]
    }
   ],
   "source": [
    "all_sv_misexp_cntrl = vrnt_id_misexp.union(vrnt_id_cntrl)\n",
    "print(f\"Total number of SVs in controls and misexpression groups: {len(all_sv_misexp_cntrl)}\")\n",
    "sv_info_misexp_cntrl_df = sv_info_df[sv_info_df.vrnt_id.isin(all_sv_misexp_cntrl)].copy()\n",
    "if sv_info_misexp_cntrl_df.shape[0] != len(all_sv_misexp_cntrl): \n",
    "    raise ValueError(\"Number of variants in SV info dataframe does not match total variants.\")\n",
    "sv_info_misexp_cntrl_df[\"group\"] = np.where(sv_info_misexp_cntrl_df.vrnt_id.isin(vrnt_id_misexp), \"Misexpression\", \"Control\")\n",
    "sv_misexp_cntrl_count_types_df = sv_info_misexp_cntrl_df.groupby([\"SVTYPE\", \"group\"], as_index=False).vrnt_id.nunique()\n",
    "sv_misexp_cntrl_count_types_df = sv_misexp_cntrl_count_types_df.rename(columns={\"vrnt_id\": \"vrnt_count\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b2861054",
   "metadata": {},
   "outputs": [],
   "source": [
    "# add labels for plot\n",
    "def label_threshold(row): \n",
    "    if row[\"group\"] == \"Control\": \n",
    "        if row[\"vrnt_count\"] > 1000: \n",
    "            return str(row[\"vrnt_count\"])\n",
    "        else:\n",
    "            return np.nan\n",
    "    else:\n",
    "        if row[\"vrnt_count\"] > 10: \n",
    "            return str(row[\"vrnt_count\"])\n",
    "        else:\n",
    "            return np.nan\n",
    "        \n",
    "sv_misexp_cntrl_count_types_df[\"label\"] = sv_misexp_cntrl_count_types_df.apply(label_threshold, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "d14b202e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# write to file \n",
    "misexp_sv_in_window_path = output_dir.joinpath(\"misexp_contrl_sv_type_count.tsv\")\n",
    "sv_misexp_cntrl_count_types_df.to_csv(misexp_sv_in_window_path, sep=\"\\t\", index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f2dc8a6",
   "metadata": {},
   "source": [
    "**Write list of selected variants and bedfiles**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "22d5467a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# write all control and misexpression-associated SVs to file \n",
    "all_sv_in_window_path = output_dir.joinpath(\"vrnt_id_in_window_cntrl_misexp_genes.txt\")\n",
    "with open(all_sv_in_window_path, \"w\") as f_out: \n",
    "    for vrnt_id in all_sv_misexp_cntrl: \n",
    "        f_out.write(f\"{vrnt_id}\\n\")\n",
    "        \n",
    "# write all misexpression-associated variants to file  \n",
    "misexp_sv_in_window_path = output_dir.joinpath(\"vrnt_id_misexp_tpm_zscore_median.txt\")\n",
    "with open(misexp_sv_in_window_path, \"w\") as f_out: \n",
    "    for vrnt_id in vrnt_id_misexp: \n",
    "        f_out.write(f\"{vrnt_id}\\n\")\n",
    "        \n",
    "# write all tested SVs to bed file (chrN format)\n",
    "vrnt_id_bed_df = sv_info_misexp_cntrl_df[[\"chr\", \"pos\", \"end\", \"vrnt_id\"]]\n",
    "vrnt_id_bed_path = output_dir.joinpath(f\"vrnt_id_in_windows_misexp_genes.bed\")\n",
    "vrnt_id_bed_df.to_csv(vrnt_id_bed_path, sep=\"\\t\", header=None, index=False)\n",
    "\n",
    "# write all tested SVs to bed file (chromosome as number format)\n",
    "vrnt_id_bed_chr_num_df = vrnt_id_bed_df.copy()\n",
    "vrnt_id_bed_chr_num_df[\"chr\"] = vrnt_id_bed_chr_num_df.chr.str.split(\"chr\").str[1]\n",
    "vrnt_id_chr_num_bed_path = output_dir.joinpath(f\"vrnt_id_in_windows_chr_num_misexp_genes.bed\")\n",
    "vrnt_id_bed_chr_num_df.to_csv(vrnt_id_chr_num_bed_path, sep=\"\\t\", header=None, index=False)\n",
    "\n",
    "# write all control and misexpression-associated SV information to file \n",
    "sv_info_misexp_cntrl_path = output_dir.joinpath(f\"vrnt_id_misexp_cntrl_info.tsv\")\n",
    "sv_info_misexp_cntrl_df.to_csv(sv_info_misexp_cntrl_path, sep=\"\\t\", index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "662203a0",
   "metadata": {},
   "source": [
    "### VEP consequence by gene "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "9e2b5142",
   "metadata": {},
   "outputs": [],
   "source": [
    "vep_msc_ranks = [\n",
    "    \"transcript_ablation\",\n",
    "    \"splice_acceptor_variant\",\n",
    "    \"splice_donor_variant\",\n",
    "    \"stop_gained\",\n",
    "    \"frameshift_variant\",\n",
    "    \"stop_lost\",\n",
    "    \"start_lost\",\n",
    "    \"transcript_amplification\",\n",
    "    \"inframe_insertion\",\n",
    "    \"inframe_deletion\",\n",
    "    \"missense_variant\",\n",
    "    \"protein_altering_variant\",\n",
    "    \"splice_donor_5th_base_variant\", \n",
    "    \"splice_region_variant\",\n",
    "    \"splice_donor_region_variant\", \n",
    "    \"splice_polypyrimidine_tract_variant\"\n",
    "    \"incomplete_terminal_codon_variant\",\n",
    "    \"start_retained_variant\",\n",
    "    \"stop_retained_variant\",\n",
    "    \"synonymous_variant\",\n",
    "    \"coding_sequence_variant\",\n",
    "    \"mature_miRNA_variant\",\n",
    "    \"5_prime_UTR_variant\",\n",
    "    \"3_prime_UTR_variant\",\n",
    "    \"non_coding_transcript_exon_variant\",\n",
    "    \"intron_variant\",\n",
    "    \"NMD_transcript_variant\",\n",
    "    \"non_coding_transcript_variant\",\n",
    "    \"coding_transcript_variant\",\n",
    "    \"upstream_gene_variant\",\n",
    "    \"downstream_gene_variant\",\n",
    "    \"TFBS_ablation\",\n",
    "    \"TFBS_amplification\",\n",
    "    \"TF_binding_site_variant\",\n",
    "    \"regulatory_region_ablation\",\n",
    "    \"regulatory_region_amplification\",\n",
    "    \"feature_elongation\",\n",
    "    \"regulatory_region_variant\",\n",
    "    \"feature_truncation\",\n",
    "    \"intergenic_variant\",\n",
    "    \"sequence_variant\",\n",
    "    \"no_predicted_effect\"\n",
    "]\n",
    "\n",
    "# from all VEP consequences that have no gene annotation\n",
    "vep_msc_not_linked_to_gene = ['TFBS_ablation', \n",
    "                              'TF_binding_site_variant',\n",
    "                              'regulatory_region_variant', \n",
    "                              'TFBS_amplification',\n",
    "                              'intergenic_variant', \n",
    "                              'regulatory_region_ablation',\n",
    "                              'regulatory_region_amplification']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "a0b44e69",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load all VEP consequences \n",
    "vep_all_df = pd.read_csv(vep_all_path, sep=\"\\t\")\n",
    "sv_gene_pairs_vep_df = pd.merge(misexp_sv_gene_pairs_df, \n",
    "                                vep_all_df[[\"vrnt_id\", \"gene_id\", \"Consequence\"]], \n",
    "                                on=[\"vrnt_id\", \"gene_id\"], \n",
    "                                how=\"left\").fillna(\"no_predicted_effect\")\n",
    "# group consequences together for each variant-gene pair\n",
    "sv_gene_pairs_vep_df[\"gene_consequence\"] = sv_gene_pairs_vep_df.groupby([\"vrnt_id\", \"gene_id\"])['Consequence'].transform(lambda x: ','.join(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "6b92b5a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove duplicates arising from variant have different consequences across gene transcripts \n",
    "sv_gene_pairs_vep_collapse_df = sv_gene_pairs_vep_df.drop(columns=[\"Consequence\"]).drop_duplicates()\n",
    "if sv_gene_pairs_vep_collapse_df.shape[0] != misexp_sv_gene_pairs_df.shape[0]: \n",
    "    raise ValueError(\"Variant gene pair number does not match number of variant gene pairs with consequence.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "284df6da",
   "metadata": {},
   "outputs": [],
   "source": [
    "# assign each variant-gene pair a unique consequence based on VEP rank  \n",
    "vrnt_gene_pair_msc_consq = []\n",
    "for index, row in sv_gene_pairs_vep_collapse_df.iterrows(): \n",
    "    gene_consequence = row[\"gene_consequence\"]\n",
    "    gene_consequence_list = gene_consequence.split(\",\")\n",
    "    for consq in vep_msc_ranks: \n",
    "        if consq in gene_consequence_list: \n",
    "            break \n",
    "    vrnt_gene_pair_msc_consq.append(consq)\n",
    "if sv_gene_pairs_vep_collapse_df.shape[0] != len(vrnt_gene_pair_msc_consq): \n",
    "    raise ValueError(\"Number of MSC per gene does not match number of variant-gene pairs.\")\n",
    "sv_gene_pairs_vep_gene_msc_df = sv_gene_pairs_vep_collapse_df.copy()\n",
    "sv_gene_pairs_vep_gene_msc_df[\"consequence\"] = vrnt_gene_pair_msc_consq\n",
    "sv_gene_pairs_vep_gene_msc_df = sv_gene_pairs_vep_gene_msc_df.drop(columns=[\"gene_consequence\"])\n",
    "# split into variants with annotated gene effect and no predicted gene effect\n",
    "vrnt_gene_effect_df = sv_gene_pairs_vep_gene_msc_df[sv_gene_pairs_vep_gene_msc_df.consequence != \"no_predicted_effect\"]\n",
    "no_predicted_effect_df = sv_gene_pairs_vep_gene_msc_df[sv_gene_pairs_vep_gene_msc_df.consequence == \"no_predicted_effect\"]\n",
    "# in cases where variant has no predicted effect add most severe consequence if regulatory \n",
    "# load VEP MSC annotations\n",
    "vep_msc_df = pd.read_csv(vep_msc_path, sep=\"\\t\").rename(columns={\"Uploaded_variation\": \"vrnt_id\", \"Consequence\": \"msc\"})\n",
    "# merge VEP MSC with no predicted effect \n",
    "no_prediced_effect_df = pd.merge(no_predicted_effect_df[[\"vrnt_id\", \"gene_id\"]], \n",
    "                                 vep_msc_df[[\"vrnt_id\", \"msc\"]], \n",
    "                                 on=\"vrnt_id\", \n",
    "                                 how=\"left\")\n",
    "# check for NaNs \n",
    "if no_prediced_effect_df.msc.isnull().values.any(): \n",
    "    raise ValueError(\"Variants missing VEP most severe consequence.\")\n",
    "no_prediced_effect_df[\"consequence\"] = np.where(no_prediced_effect_df.msc.isin(vep_msc_not_linked_to_gene), \n",
    "                                                no_prediced_effect_df.msc, \n",
    "                                                \"no_predicted_effect\")\n",
    "no_prediced_effect_trunc_df = no_prediced_effect_df[[\"vrnt_id\", \"gene_id\", \"consequence\"]]\n",
    "# combine variants with annotated gene effect and variants with no predicted effect with updated regulatory consequence\n",
    "vrnt_gene_pair_consq_msc_added_df = pd.concat([vrnt_gene_effect_df, no_prediced_effect_trunc_df])\n",
    "if vrnt_gene_pair_consq_msc_added_df[[\"vrnt_id\", \"gene_id\"]].drop_duplicates().shape[0] != misexp_sv_gene_pairs_df.shape[0]: \n",
    "    raise ValueError(\"Number of gene pairs in expression input does not match number of gene pairs with annotated consequence.\")\n",
    "# add SV type\n",
    "vrnt_gene_pair_consq_msc_added_df = pd.merge(vrnt_gene_pair_consq_msc_added_df, \n",
    "                                             sv_info_df[[\"vrnt_id\", \"SVTYPE\"]], \n",
    "                                             on=\"vrnt_id\", \n",
    "                                             how=\"inner\")\n",
    "# write to file \n",
    "vrnt_id_gene_consq_msc_path = output_dir.joinpath(\"misexp_vrnt_gene_msc_consq.tsv\")\n",
    "vrnt_gene_pair_consq_msc_added_df.to_csv(vrnt_id_gene_consq_msc_path, sep=\"\\t\", index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "b1cb8662",
   "metadata": {},
   "outputs": [],
   "source": [
    "# consequence names \n",
    "consequence_names_dict = {\n",
    "    '3_prime_UTR_variant': \"3' UTR variant\",\n",
    "    '5_prime_UTR_variant': \"5' UTR variant\",\n",
    "    'TFBS_ablation': \"TFBS ablation\",\n",
    "    'TF_binding_site_variant': \"TF binding site variant\",\n",
    "    'coding_sequence_variant': \"Coding sequence variant\",\n",
    "    'downstream_gene_variant': \"Downstream gene variant\",\n",
    "    'frameshift_variant': \"Frameshift variant\",\n",
    "    'inframe_deletion': \"Inframe deletion\",\n",
    "    'intergenic_variant': \"Intergenic variant\",\n",
    "    'intron_variant': \"Intron variant\",\n",
    "    'no_predicted_effect': \"No predicted effect on gene\",\n",
    "    'non_coding_transcript_exon_variant': \"Non-coding transcript exon variant\",\n",
    "    'regulatory_region_ablation': \"Regulatory region ablation\",\n",
    "    'regulatory_region_variant': \"Regulatory region variant\",\n",
    "    'stop_lost': \"Stop lost\",\n",
    "    'transcript_ablation': \"Transcript ablation\",\n",
    "    'upstream_gene_variant': \"Upstream gene variant\",\n",
    "    'TFBS_amplification': \"TFBS amplification\",\n",
    "    'regulatory_region_amplification': \"Regulatory region amplification\",\n",
    "    'transcript_amplification': \"Transcript amplification\"\n",
    "}\n",
    "\n",
    "# consequence groups \n",
    "conseq_group_dict = {\n",
    "    '3_prime_UTR_variant': \"3' UTR variant\",\n",
    "    '5_prime_UTR_variant': \"5' UTR variant\",\n",
    "    'TFBS_ablation': \"Regulatory\",\n",
    "    'TF_binding_site_variant': \"Regulatory\",\n",
    "    'coding_sequence_variant': \"Coding\",\n",
    "    'downstream_gene_variant': \"Downstream (5 kb)\",\n",
    "    'frameshift_variant': \"Coding\",\n",
    "    'inframe_deletion': \"Coding\",\n",
    "    'intergenic_variant': \"Non-coding\",\n",
    "    'intron_variant': \"Non-coding\",\n",
    "    'no_predicted_effect': \"No predicted effect\",\n",
    "    'non_coding_transcript_exon_variant': \"Non-coding\",\n",
    "    'regulatory_region_ablation': \"Regulatory\",\n",
    "    'regulatory_region_variant': \"Regulatory\",\n",
    "    'stop_lost': \"Coding\",\n",
    "    'transcript_ablation': \"Transcript ablation\",\n",
    "    'upstream_gene_variant': \"Upstream (5 kb)\",\n",
    "    'TFBS_amplification': \"Regulatory\",\n",
    "    'regulatory_region_amplification': \"Regulatory\",\n",
    "    'transcript_amplification': \"Transcript amplification\"    \n",
    "}\n",
    "\n",
    "# count variant consequences\n",
    "vrnt_id_gene_conseq_count_df = vrnt_gene_pair_consq_msc_added_df.groupby([\"SVTYPE\", \"consequence\"], as_index=False).vrnt_id.count()\n",
    "vrnt_id_gene_conseq_count_df[\"consq_name\"] = vrnt_id_gene_conseq_count_df.consequence.replace(consequence_names_dict)\n",
    "vrnt_id_gene_conseq_count_df[\"consq_group\"] = vrnt_id_gene_conseq_count_df.consequence.replace(conseq_group_dict)\n",
    "# check for NaNs \n",
    "if not vrnt_id_gene_conseq_count_df[vrnt_id_gene_conseq_count_df[\"consq_name\"].isna()].empty: \n",
    "    raise ValueError(\"Contains NaNs - missing consequence names.\")\n",
    "if not vrnt_id_gene_conseq_count_df[vrnt_id_gene_conseq_count_df[\"consq_group\"].isna()].empty: \n",
    "    raise ValueError(\"Contains NaNs - missing consequence groups.\")\n",
    "vrnt_id_gene_conseq_count_df = vrnt_id_gene_conseq_count_df.rename(columns={\"vrnt_id\": \"vrnt_gene_count\"})\n",
    "\n",
    "# calculate proportion of each variant consequence relative to total in SVTYPE group \n",
    "count_vrnt_in_group_df = vrnt_id_gene_conseq_count_df.groupby([\"SVTYPE\"], as_index=False).vrnt_gene_count.sum()\n",
    "count_vrnt_in_group_df = count_vrnt_in_group_df.rename(columns={\"vrnt_gene_count\":\"vrnt_gene_in_svtype_group\"})\n",
    "vrnt_id_gene_conseq_count_proport_df = pd.merge(vrnt_id_gene_conseq_count_df, \n",
    "                                        count_vrnt_in_group_df, \n",
    "                                        on=[\"SVTYPE\"], \n",
    "                                        how=\"inner\")\n",
    "vrnt_id_gene_conseq_count_proport_df[\"proport\"] = vrnt_id_gene_conseq_count_proport_df.vrnt_gene_count / vrnt_id_gene_conseq_count_proport_df.vrnt_gene_in_svtype_group"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "d5864004",
   "metadata": {},
   "outputs": [],
   "source": [
    "# write to file \n",
    "vrnt_id_gene_conseq_count_path = output_dir.joinpath(\"vrnt_consq_count.tsv\")\n",
    "vrnt_id_gene_conseq_count_proport_df.to_csv(vrnt_id_gene_conseq_count_path, sep=\"\\t\", index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "60032ab8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# counts of only deletions and duplications \n",
    "del_dup_gene_conseq_count_proport_df = vrnt_id_gene_conseq_count_proport_df[vrnt_id_gene_conseq_count_proport_df.SVTYPE.isin([\"DEL\", \"DUP\"])]\n",
    "del_dup_gene_conseq_count_proport_path = output_dir.joinpath(\"del_dup_consq_count.tsv\")\n",
    "del_dup_gene_conseq_count_proport_df.to_csv(del_dup_gene_conseq_count_proport_path, sep=\"\\t\", index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "b67e0de3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# count consequences by group \n",
    "del_dup_gene_conseq_group_count_proport_df = del_dup_gene_conseq_count_proport_df.groupby([\"SVTYPE\", \"consq_group\"], as_index=False).vrnt_gene_count.sum()\n",
    "del_dup_gene_conseq_group_count_proport_path = output_dir.joinpath(\"del_dup_consq_group_count.tsv\")\n",
    "del_dup_gene_conseq_group_count_proport_df.to_csv(del_dup_gene_conseq_group_count_proport_path, sep=\"\\t\", index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "398fd2fb",
   "metadata": {},
   "source": [
    "### Positions of misexpression-associated SVs relative to misexpressed gene "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "6e227a1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# add SV info \n",
    "vrnt_gene_pair_consq_msc_info_df = pd.merge(vrnt_gene_pair_consq_msc_added_df,\n",
    "                                            sv_info_df[[\"vrnt_id\", \"chr\", \"pos\", \"end\", \"SVLEN\"]], \n",
    "                                            on=\"vrnt_id\", \n",
    "                                            how=\"inner\"\n",
    "                                           )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "d1ab6d1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# add gene position \n",
    "gene_ids_misexp_cntrl = vrnt_gene_pair_consq_msc_info_df.gene_id.unique()\n",
    "overlap_gene_coord_dict, count = {}, 0\n",
    "for gtf in pysam.TabixFile(str(gencode_path)).fetch(parser = pysam.asGTF()): \n",
    "    if gtf.feature == \"gene\" and gtf.gene_id.split(\".\")[0] in gene_ids_misexp_cntrl: \n",
    "        overlap_gene_coord_dict[count] = [gtf.gene_id.split(\".\")[0], gtf.contig, gtf.start, gtf.end, gtf.strand]\n",
    "        count += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "feae2012",
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_coord_colnames = [\"gene_id\", \"gene_chrom\", \"gene_start\", \"gene_end\", \"gene_strand\"]\n",
    "gene_coord_df = pd.DataFrame.from_dict(overlap_gene_coord_dict, orient=\"index\", columns=gene_coord_colnames)\n",
    "vrnt_gene_pair_consq_msc_info_gene_df = pd.merge(vrnt_gene_pair_consq_msc_info_df, \n",
    "                                                 gene_coord_df, \n",
    "                                                 on=[\"gene_id\"], \n",
    "                                                 how=\"inner\")\n",
    "# check chromosomes match\n",
    "if not vrnt_gene_pair_consq_msc_info_gene_df[vrnt_gene_pair_consq_msc_info_gene_df.chr != vrnt_gene_pair_consq_msc_info_gene_df.gene_chrom].empty:\n",
    "    raise ValueError(\"Chromosomes do not match\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "e723a39d",
   "metadata": {},
   "outputs": [],
   "source": [
    "### check whether misexpressed gene is overlapped by variant \n",
    "vrnt_gene_pairs_df = vrnt_gene_pair_consq_msc_info_gene_df[[\"vrnt_id\", \"gene_id\", \"chr\", \"pos\", \"end\", \"SVTYPE\", \"consequence\", \"gene_start\", \"gene_end\", \"gene_strand\"]].drop_duplicates()\n",
    "vrnt_id_pos_to_gene, count = {}, 0\n",
    "for index, row in vrnt_gene_pairs_df.iterrows(): \n",
    "    vrnt_id, gene_id, msc = row[\"vrnt_id\"], row[\"gene_id\"], row[\"consequence\"]\n",
    "    gene_strand, sv_type = row[\"gene_strand\"], row[\"SVTYPE\"]\n",
    "    sv_start, sv_end =  row[\"pos\"], row[\"end\"]\n",
    "    gene_chrom, gene_start, gene_end = row[\"chr\"], row[\"gene_start\"], row[\"gene_end\"]\n",
    "    if gene_strand == \"+\":\n",
    "        # upstream \n",
    "        if sv_start < gene_start and sv_end < gene_start: \n",
    "            position = \"Upstream\"\n",
    "        # downstream\n",
    "        elif sv_start > gene_end and sv_end > gene_end: \n",
    "            position = \"Downstream\"\n",
    "        # internal \n",
    "        elif sv_start > gene_start and gene_end > sv_end:\n",
    "            position = \"Internal\"\n",
    "        # entire gene \n",
    "        elif sv_start < gene_start and gene_end < sv_end: \n",
    "            position = \"Entire gene\"\n",
    "        # 5 prime overlap\n",
    "        elif sv_start <= gene_start and sv_end > gene_start and gene_end > sv_end:  \n",
    "            position = \"Partial overlap 5' end\"\n",
    "        # 3 prime overlap \n",
    "        elif sv_start < gene_end and sv_start > gene_start and sv_end >= gene_end: \n",
    "            position = \"Partial overlap 3' end\"\n",
    "        else: \n",
    "            raise ValueError(f\"{vrnt_id} {gene_id} no position found.\")\n",
    "    if gene_strand == \"-\": \n",
    "        # upstream \n",
    "        if sv_start > gene_end and sv_end > gene_end: \n",
    "            position = \"Upstream\"\n",
    "        # downstream \n",
    "        elif sv_start < gene_start and sv_end < gene_start: \n",
    "            position = \"Downstream\"\n",
    "        # internal \n",
    "        elif sv_start > gene_start and gene_end > sv_end:\n",
    "            position = \"Internal\"\n",
    "        # entire gene \n",
    "        elif sv_start <= gene_start and gene_end <= sv_end: \n",
    "            position = \"Entire gene\"\n",
    "        # 5 prime overlap \n",
    "        elif sv_start > gene_start and gene_end > sv_start and gene_end <= sv_end: \n",
    "            position = \"Partial overlap 5' end\"\n",
    "        # 3 prime overlap\n",
    "        elif sv_start <= gene_start and gene_end > sv_end and gene_start < sv_end:  \n",
    "            position = \"Partial overlap 3' end\"\n",
    "        else: \n",
    "            raise ValueError(f\"{vrnt_id} {gene_id} no position found.\")\n",
    "    vrnt_id_pos_to_gene[count] = [gene_chrom, vrnt_id, sv_start, sv_end, gene_id, \n",
    "                                         gene_start, gene_end, gene_strand, position, msc, sv_type]\n",
    "    count += 1\n",
    "    \n",
    "vrnt_id_pos_cols = [\"chrom\", \"vrnt_id\", \"sv_start\", \"sv_end\", \"gene_id\", \"gene_start\", \n",
    "                    \"gene_end\", \"gene_strand\", \"position\", \"consequence\", \"SVTYPE\"]\n",
    "vrnt_id_pos_to_gene_df = pd.DataFrame.from_dict(vrnt_id_pos_to_gene, \n",
    "                                                       orient=\"index\", \n",
    "                                                       columns=vrnt_id_pos_cols)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "e619ba54",
   "metadata": {},
   "outputs": [],
   "source": [
    "vrnt_id_pos_to_gene_path = output_dir.joinpath(\"misexp_vrnt_gene_position.tsv\")\n",
    "vrnt_id_pos_to_gene_df.to_csv(vrnt_id_pos_to_gene_path, sep=\"\\t\", index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "57fb5a83",
   "metadata": {},
   "outputs": [],
   "source": [
    "# count positions of misexpression-associated variants  \n",
    "vrnt_id_pos_to_gene_count_df = vrnt_id_pos_to_gene_df.groupby([\"SVTYPE\", \"position\"], as_index=False).vrnt_id.count()\n",
    "vrnt_id_pos_to_gene_count_df = vrnt_id_pos_to_gene_count_df.rename(columns={\"vrnt_id\": \"vrnt_count\"})\n",
    "# add proportion\n",
    "vrnt_id_pos_total_df = vrnt_id_pos_to_gene_count_df.groupby([\"SVTYPE\"], as_index=False).vrnt_count.sum()\n",
    "vrnt_id_pos_total_df = vrnt_id_pos_total_df.rename(columns={\"vrnt_count\": \"total_in_group\"})\n",
    "vrnt_id_pos_to_gene_proport_df = pd.merge(vrnt_id_pos_to_gene_count_df, \n",
    "                                          vrnt_id_pos_total_df, \n",
    "                                          on=[\"SVTYPE\"], \n",
    "                                          how=\"inner\")\n",
    "vrnt_id_pos_to_gene_proport_df[\"proport\"] = vrnt_id_pos_to_gene_proport_df.vrnt_count / vrnt_id_pos_to_gene_proport_df.total_in_group\n",
    "# write to file \n",
    "vrnt_id_pos_to_gene_proport_path = output_dir.joinpath(\"vrnt_id_pos_to_gene_proport.tsv\")\n",
    "vrnt_id_pos_to_gene_proport_df.to_csv(vrnt_id_pos_to_gene_proport_path, sep=\"\\t\", index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "c0dcf9fb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Percentage of DUPs overlap the misexpressed gene: 72.22222222222223\n",
      "Percentage of DELs that do overlap the misexpressed gene: 7.526881720430108\n"
     ]
    }
   ],
   "source": [
    "# proportion of DUPs overlapping the misexpressed gene entirely or partially \n",
    "partial_entire_overlap = [\"Entire gene\", \"Partial overlap 3' end\", \"Partial overlap 5' end\" , \"internal\"]\n",
    "proport_dup_gene = vrnt_id_pos_to_gene_proport_df[(vrnt_id_pos_to_gene_proport_df.position.isin(partial_entire_overlap)) & \n",
    "                                                  (vrnt_id_pos_to_gene_proport_df.SVTYPE == \"DUP\")\n",
    "                                                 ].proport.sum()\n",
    "print(f\"Percentage of DUPs overlap the misexpressed gene: {proport_dup_gene*100}\")\n",
    "proport_del_not_gene = vrnt_id_pos_to_gene_proport_df[(vrnt_id_pos_to_gene_proport_df.position.isin(partial_entire_overlap)) & \n",
    "                                                  (vrnt_id_pos_to_gene_proport_df.SVTYPE == \"DEL\")\n",
    "                                                 ].proport.sum()\n",
    "print(f\"Percentage of DELs that do overlap the misexpressed gene: {proport_del_not_gene*100}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "f26d2c93",
   "metadata": {},
   "outputs": [],
   "source": [
    "# write information to file \n",
    "vrnt_gene_pair_consq_msc_info_gene_pos_df = pd.merge(vrnt_gene_pair_consq_msc_info_gene_df, \n",
    "                                                 vrnt_id_pos_to_gene_df[[\"vrnt_id\", \"gene_id\", \"position\"]], \n",
    "                                                 on=[\"vrnt_id\", \"gene_id\"], \n",
    "                                                 how=\"inner\")\n",
    "vrnt_gene_pair_consq_msc_info_gene_pos_path = output_dir.joinpath(\"vrnt_gene_pair_misexp_cntrl_info.tsv\")\n",
    "vrnt_gene_pair_consq_msc_info_gene_pos_df.to_csv(vrnt_gene_pair_consq_msc_info_gene_pos_path, sep=\"\\t\", index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "62323fa3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# write deletions and duplications only to file \n",
    "dup_del_gene_pair_consq_msc_info_gene_pos_df = vrnt_gene_pair_consq_msc_info_gene_pos_df[vrnt_gene_pair_consq_msc_info_gene_pos_df.SVTYPE.isin([\"DEL\", \"DUP\"])]\n",
    "dup_del_gene_pair_consq_msc_info_gene_pos_path = output_dir.joinpath(\"dup_del_gene_pair_misexp_cntrl_info.tsv\")\n",
    "dup_del_gene_pair_consq_msc_info_gene_pos_df.to_csv(dup_del_gene_pair_consq_msc_info_gene_pos_path, sep=\"\\t\", index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "db6b13ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "### comparison of position and VEP prediction "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "569c8fe8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SVTYPE</th>\n",
       "      <th>consequence</th>\n",
       "      <th>position</th>\n",
       "      <th>vrnt_id</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>DEL</td>\n",
       "      <td>TFBS_ablation</td>\n",
       "      <td>Upstream</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>DEL</td>\n",
       "      <td>coding_sequence_variant</td>\n",
       "      <td>Partial overlap 5' end</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>DEL</td>\n",
       "      <td>downstream_gene_variant</td>\n",
       "      <td>Downstream</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>DEL</td>\n",
       "      <td>intergenic_variant</td>\n",
       "      <td>Upstream</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>DEL</td>\n",
       "      <td>intron_variant</td>\n",
       "      <td>Internal</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>DEL</td>\n",
       "      <td>no_predicted_effect</td>\n",
       "      <td>Downstream</td>\n",
       "      <td>21</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>DEL</td>\n",
       "      <td>no_predicted_effect</td>\n",
       "      <td>Upstream</td>\n",
       "      <td>35</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>DEL</td>\n",
       "      <td>non_coding_transcript_exon_variant</td>\n",
       "      <td>Partial overlap 5' end</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>DEL</td>\n",
       "      <td>regulatory_region_ablation</td>\n",
       "      <td>Upstream</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>DEL</td>\n",
       "      <td>regulatory_region_variant</td>\n",
       "      <td>Downstream</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>DEL</td>\n",
       "      <td>upstream_gene_variant</td>\n",
       "      <td>Upstream</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>DUP</td>\n",
       "      <td>5_prime_UTR_variant</td>\n",
       "      <td>Partial overlap 5' end</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>DUP</td>\n",
       "      <td>coding_sequence_variant</td>\n",
       "      <td>Partial overlap 3' end</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>DUP</td>\n",
       "      <td>no_predicted_effect</td>\n",
       "      <td>Downstream</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>DUP</td>\n",
       "      <td>no_predicted_effect</td>\n",
       "      <td>Upstream</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>DUP</td>\n",
       "      <td>non_coding_transcript_exon_variant</td>\n",
       "      <td>Entire gene</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>DUP</td>\n",
       "      <td>non_coding_transcript_exon_variant</td>\n",
       "      <td>Partial overlap 3' end</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>DUP</td>\n",
       "      <td>transcript_amplification</td>\n",
       "      <td>Entire gene</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>DUP</td>\n",
       "      <td>transcript_amplification</td>\n",
       "      <td>Partial overlap 3' end</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>INV</td>\n",
       "      <td>coding_sequence_variant</td>\n",
       "      <td>Partial overlap 5' end</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>INV</td>\n",
       "      <td>no_predicted_effect</td>\n",
       "      <td>Upstream</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   SVTYPE                         consequence                position  vrnt_id\n",
       "0     DEL                       TFBS_ablation                Upstream        2\n",
       "1     DEL             coding_sequence_variant  Partial overlap 5' end        3\n",
       "2     DEL             downstream_gene_variant              Downstream        1\n",
       "3     DEL                  intergenic_variant                Upstream        9\n",
       "4     DEL                      intron_variant                Internal        5\n",
       "5     DEL                 no_predicted_effect              Downstream       21\n",
       "6     DEL                 no_predicted_effect                Upstream       35\n",
       "7     DEL  non_coding_transcript_exon_variant  Partial overlap 5' end        4\n",
       "8     DEL          regulatory_region_ablation                Upstream        3\n",
       "9     DEL           regulatory_region_variant              Downstream        3\n",
       "10    DEL               upstream_gene_variant                Upstream        7\n",
       "11    DUP                 5_prime_UTR_variant  Partial overlap 5' end        1\n",
       "12    DUP             coding_sequence_variant  Partial overlap 3' end        3\n",
       "13    DUP                 no_predicted_effect              Downstream        2\n",
       "14    DUP                 no_predicted_effect                Upstream        3\n",
       "15    DUP  non_coding_transcript_exon_variant             Entire gene        3\n",
       "16    DUP  non_coding_transcript_exon_variant  Partial overlap 3' end        1\n",
       "17    DUP            transcript_amplification             Entire gene        4\n",
       "18    DUP            transcript_amplification  Partial overlap 3' end        1\n",
       "19    INV             coding_sequence_variant  Partial overlap 5' end        1\n",
       "20    INV                 no_predicted_effect                Upstream        1"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "misexp_vrnt_gene_pair_consq_msc_info_gene_pos_df = vrnt_gene_pair_consq_msc_info_gene_pos_df\n",
    "misexp_vrnt_gene_pair_consq_msc_info_gene_pos_df.groupby([\"SVTYPE\", \"consequence\", \"position\"], as_index=False).vrnt_id.count()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f970b9b",
   "metadata": {},
   "source": [
    "### Comparison of SV lengths (all SVs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "162bbb2e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DEL one-side Mann-whitney length difference p-value: 6.029879821698881e-07\n",
      "Control DEL mean length: 9745.010125307675\n",
      "Misexpression DEL mean length: 24515.011494252874\n",
      "Average fold-change in DEL length: 2.5156476164747823\n",
      "DUP one-side Mann-whitney length difference p-value: 3.9317030193917875e-06\n",
      "Control DUP mean length: 58723.69349112426\n",
      "Misexpression DUP mean length: 199619.3125\n",
      "Average fold-change in DUP length: 3.399297636654467\n"
     ]
    }
   ],
   "source": [
    "# one-sided Mann-Whitney U-test comparison of control and misexpression length\n",
    "for sv_type in [\"DEL\", \"DUP\"]: \n",
    "    sv_type_info_misexp_cntrl_df = sv_info_misexp_cntrl_df[sv_info_misexp_cntrl_df.SVTYPE == sv_type]\n",
    "    length_controls = sv_type_info_misexp_cntrl_df[sv_type_info_misexp_cntrl_df.group == \"Control\"].SVLEN.to_numpy()\n",
    "    length_misexp = sv_type_info_misexp_cntrl_df[sv_type_info_misexp_cntrl_df.group == \"Misexpression\"].SVLEN.to_numpy()\n",
    "    # ‘less’: the distribution underlying x is stochastically less than the distribution \n",
    "    # underlying y, i.e. F(u) > G(u) for all u.\n",
    "    stat, pval = mannwhitneyu(length_controls, length_misexp, alternative=\"less\")\n",
    "    print(f\"{sv_type} one-side Mann-whitney length difference p-value: {pval}\")\n",
    "    # fold change in length \n",
    "    mean_length_controls = np.average(length_controls)\n",
    "    print(f\"Control {sv_type} mean length: {mean_length_controls}\")\n",
    "    mean_length_misexp = np.average(length_misexp)\n",
    "    print(f\"Misexpression {sv_type} mean length: {mean_length_misexp}\")\n",
    "    print(f\"Average fold-change in {sv_type} length: {mean_length_misexp/mean_length_controls}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "2964eaaf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DEL one-side Mann-whitney length difference p-value: 0.00023181779155355376\n",
      "Control DEL mean length: 14621.333463110765\n",
      "Misexpression DEL mean length: 41628.117647058825\n",
      "Fold-change in DEL length: 2.8470807913714204\n",
      "DUP one-side Mann-whitney length difference p-value: 0.0022695242582722016\n",
      "Control DUP mean length: 72434.5397225725\n",
      "Misexpression DUP mean length: 144444.36363636365\n",
      "Fold-change in DUP length: 1.9941365568082845\n"
     ]
    }
   ],
   "source": [
    "# one-sided Mann-Whitney U-test (restricted to singletons)\n",
    "for sv_type in [\"DEL\", \"DUP\"]: \n",
    "    sv_type_info_misexp_cntrl_df = sv_info_misexp_cntrl_df[(sv_info_misexp_cntrl_df.SVTYPE == sv_type) & \n",
    "                                                           (sv_info_misexp_cntrl_df.AC == 1)\n",
    "                                                          ]\n",
    "    length_controls = sv_type_info_misexp_cntrl_df[sv_type_info_misexp_cntrl_df.group == \"Control\"].SVLEN.to_numpy()\n",
    "    length_misexp = sv_type_info_misexp_cntrl_df[sv_type_info_misexp_cntrl_df.group == \"Misexpression\"].SVLEN.to_numpy()\n",
    "    # ‘less’: the distribution underlying x is stochastically less than the distribution \n",
    "    # underlying y, i.e. F(u) > G(u) for all u.\n",
    "    stat, pval = mannwhitneyu(length_controls, length_misexp, alternative=\"less\")\n",
    "    print(f\"{sv_type} one-side Mann-whitney length difference p-value: {pval}\")\n",
    "    # fold change in length \n",
    "    mean_length_controls = np.average(length_controls)\n",
    "    print(f\"Control {sv_type} mean length: {mean_length_controls}\")\n",
    "    mean_length_misexp = np.average(length_misexp)\n",
    "    print(f\"Misexpression {sv_type} mean length: {mean_length_misexp}\")\n",
    "    print(f\"Fold-change in {sv_type} length: {mean_length_misexp/mean_length_controls}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bc3af28",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
